
alpha = -180:5:180;
beta = -90:5:90;

alpha_r = deg2rad(alpha);
beta_r = deg2rad(beta);

[alpha_mesh, beta_mesh] = meshgrid(alpha_r, beta_r);


beta_f=0;
alpha_f=30;

beta_f1 = deg2rad(beta_f);
alpha_f1 = deg2rad(alpha_f);

R=[1/4 -1/4 0; -1/4 -1/4 0; 1/4 -1/4 0; 1/4 1/4 0];
K=[0 cos(alpha_f) sin(alpha_f)];
dPhi=2*pi*K*R';
B = exp(-1i*dPhi);
for i = 1: size(alpha_mesh,1);
    for j = 1: size(alpha_mesh,2);    
alpha = alpha_mesh(i,j);
beta = beta_mesh(i,j);
k = [ cos(alpha)*sin(beta) sin(alpha)*sin(beta) cos(beta)];
B1=exp(+1*2*pi*k*R')';
l(i,j)=abs(sum(B*B1))^2;
end;
end;
    
[x, y, z] = sph2cart(beta_mesh, alpha_mesh, l);

figure(1);
surf( x, y, z)
xlabel('x'); ylabel('y'); zlabel('z')

minc = min( [min(min(x)) min(min(y)) min(min(z))] );
maxc = max( [max(max(x)) max(max(y)) max(max(z))] );
xlim([minc maxc]); ylim([minc maxc]); zlim([minc maxc]);

% surf( alpha_mesh, beta_mesh, R)